library(foreign)
library(dplyr)
library(extrafont)

data = data.frame(read.dta("C:/Users/Michael/Documents/fop/psp_v/revisions/2ndRevisions/datacode/FINAL/s1/fopdata.dta"))

state_df = data %>%
  group_by(state) %>%
  summarise(TotalVotes = sum(total_2016),TotalTrumpVotes = sum(gop_2016), TotalClintonVotes = sum(dem_2016))
state_df = data.frame(state_df)
state_df$ClintonVS = state_df$TotalClintonVotes/state_df$TotalVotes
state_df$TrumpVS = state_df$TotalTrumpVotes/state_df$TotalVotes

# given beta, generates cf elections 
getstates = function(data, beta) {
  holder = data
  holder$cf = data$pct_trump - (beta/100*data$log_lodge_density)
  holder$CF_Vote = holder$cf * data$total_2016
  
  hstate_df = holder %>% group_by(state) %>% summarise(TrumpCF = sum(CF_Vote), TotalClintonVotes = sum(dem_2016),
                                                       TotalVotes = sum(total_2016))
  hstate_df$CF_TVS = hstate_df$TrumpCF/hstate_df$TotalVotes
  hstate_df$CF_Margin = hstate_df$CF_TVS-(hstate_df$TotalClintonVotes/hstate_df$TotalVotes)
  return(hstate_df[,c("CF_TVS", "CF_Margin")])
}

set.seed(02138)
NREP = 10000
dist = rnorm(NREP, mean = .8898729, sd=  .2259673)

my_results <- list()
idx = 1
for(beta in dist) {
  my_results[[idx]] = getstates(data, beta)
  idx = idx+1
}
my_results = lapply(my_results, data.frame)

lapply(lapply(my_results, "[[", 1), "[[", 1)

stateplot = function(ST, state_df, xa=0.478, ya=0.4595, za=0.461) {
  if(state_df[which(state_df$state==ST), "ClintonVS"] < state_df[which(state_df$state==ST), "TrumpVS"]) {
    a = plot(density(unlist(lapply(lapply(my_results, "[[", 1), "[[", which(state_df$state==ST)))), xlim=c(state_df[which(state_df$state==ST), "ClintonVS"]-0.01, state_df[which(state_df$state==ST), "TrumpVS"]+0.01),
             lty='solid', col="pink", lwd=1, cex.main = 2, xaxt="n", xlab="Vote Share" ,main = paste("Counterfactual Election in", ST))
    a = axis(1, at=pretty(c(state_df[which(state_df$state==ST), "ClintonVS"]-0.01, state_df[which(state_df$state==ST), "TrumpVS"]+0.01)), lab=paste0(pretty(c(state_df[which(state_df$state==ST), "ClintonVS"]-0.01, state_df[which(state_df$state==ST), "TrumpVS"]+0.01)) * 100, "%"),
             cex.axis = 1.4, las=TRUE)
    a = a+ polygon(density(unlist(lapply(lapply(my_results, "[[", 1), "[[", which(state_df$state==ST)))), col=adjustcolor("red",alpha.f=0.3), border=NA, lwd=3)
    a = a + abline(v=state_df[which(state_df$state==ST), "ClintonVS"], col="blue", lty="dotted", lwd=3)
    a = a + abline(v=state_df[which(state_df$state==ST), "TrumpVS"], col="red", lty="dotted", lwd=3)
    pr = mean(unlist(lapply(lapply(my_results, "[[", 1), "[[", which(state_df$state==ST)))<state_df[which(state_df$state==ST), "ClintonVS"])
    print(pr)
    legend("right", legend = c("Trump", "Clinton"),
                   bty="n",col=c("red", "blue"), lty=c(2,2), cex=1.4)
    legend(xa, 200, legend = c("No FOP"),
           bty="n", fill=adjustcolor("red",alpha.f=0.3), cex=1.4)
    if(pr > 0) {
      a = a + text(ya, 180 , paste("Pr(Pivotal)", " = ", round(pr,3), sep = ""), 
           adj = 0, cex = 1.6)
    }
    a = a + text(za, 100 , paste("(",round(sort(ci_es(ST, state_df,
                      ci=0.99))[1], 2), "%, ",round(sort(ci_es(ST, state_df,
                      ci=0.99))[2], 2), "%)", sep = ""), adj = 0, cex = 1.3)
    a = a + text(za, 125 , paste("Effect, 99% CI", sep = ""), adj = 0, cex = 1.3)

    return(a)
  }
  else {
    a = plot(density(unlist(lapply(lapply(my_results, "[[", 1), "[[", which(state_df$state==ST)))), xlim=c(state_df[which(state_df$state==ST), "TrumpVS"]-0.01, state_df[which(state_df$state==ST), "ClintonVS"]+0.01),
             lty='solid', col="pink", lwd=1, xaxt="n", xlab="Vote Share", cex.main = 2, main = paste("Counterfactual Election in", ST))
    a = a + abline(v=state_df[which(state_df$state==ST), "ClintonVS"], col="blue", lty="dotted", lwd=3)
    a = axis(1, at=pretty(c(state_df[which(state_df$state==ST), "TrumpVS"]-0.01, state_df[which(state_df$state==ST), "ClintonVS"]+0.01)), lab=paste0(pretty(c(state_df[which(state_df$state==ST), "ClintonVS"]+0.01, state_df[which(state_df$state==ST), "TrumpVS"]-0.01)) * 100, "%"),
             cex.axis = 1.4,las=TRUE)
    a = a + polygon(density(unlist(lapply(lapply(my_results, "[[", 1), "[[", which(state_df$state==ST)))), col=adjustcolor("red",alpha.f=0.3), border=NA, lwd=3)
    a = a + abline(v=state_df[which(state_df$state==ST), "TrumpVS"], col="red", lty="dotted", lwd=3)
    legend("topright", legend = c("Trump", "Clinton"),
           bty="n",col=c("red", "blue"), lty=c(2,2), cex=1.4)
    legend("right", legend = c("No FOP"),
           bty="n", fill=adjustcolor("red",alpha.f=0.3), cex=1.4)
    return(a)
  }
  a
}

ci_es = function(ST, state_df, ci=0.99) {
  return((state_df[which(state_df$state==ST), "TrumpVS"]-quantile(unlist(lapply(lapply(my_results,
                                                                                       "[[", 1), "[[", which(state_df$state==ST))), c(0.5-0.5*ci, 0.5+0.5*ci)))*100)
}

setwd("C:/Users/Michael/Documents/fop/psp_v/revisions/2ndRevisions/datacode/")

#stateplot("NH", state_df)
png('mi_cf.png', units="in", width=8, height=8, res=300)
stateplot("MI", state_df)
dev.off()

png('pa_cf.png', units="in", width=8, height=8, res=300)
stateplot("PA", state_df, 0.485, 0.465, 0.467)
dev.off()

png('nh_cf.png', units="in", width=8, height=8, res=300)
stateplot("NH", state_df)
dev.off()

stateplot("FL", state_df)
stateplot("NC", state_df)
stateplot("VA", state_df)
stateplot("WI", state_df, 0.475, 0.458, 0.457)

plots = list()
plots = lapply(state_df$state, stateplot, state_df)


med_es = function(ST, state_df) {
  return((state_df[which(state_df$state==ST), "TrumpVS"]-quantile(unlist(lapply(lapply(my_results,
      "[[", 1), "[[", which(state_df$state==ST))), c(0.5)))*100)
}

mean_es = function(ST, state_df) {
  return((state_df[which(state_df$state==ST), "TrumpVS"]-mean(unlist(lapply(lapply(my_results, "[[", 1), "[[", which(state_df$state==ST)))))*100)
}

margin = function(ST, state_df) {
  return((state_df[which(state_df$state==ST), "TrumpVS"]-state_df[which(state_df$state==ST), "ClintonVS"])*100)
}

pivotality = function(ST, state_df) {
  abs(mean_es(ST, state_df)/margin(ST, state_df))
}

state_df["Pivotality"] = unlist(lapply(state_df$state, pivotality, state_df))

